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Calculation of the nonrelativistic Bethe logarithm in the velocity gauge 

V.I. Korobov 
Joint Institute for Nuclear Research, 141980, Dubna, Russia 

We consider a general procedure to evaluate the Bethe logarithm for a general few-body atomic 
or molecular system. As benchmarks we use calculation for the ground states of a helium atom 
and Hj molecular ion. The obtained values are: /3h c = 4.37016022306(2) for the helium atom and 
/3 H + = 3.012230335(1) for the Hj. Both results substantially improve the best known values for 
these quantities. 



O ■ I. INTRODUCTION 



It is known that in the leading order radiative contribution to the binding energy of atomic or molecular system the 
most complicate for numerical evaluation quantity is the Bethe logarithm [l| . One of the first most accurate results 
for the helium ground state was obtained by C. Schwartz Q in 1961 and remained the best one over 30 years! In 1999, 
CN . Goldman and Drake [3| suggested a new way to evaluate the Bethe logarithm, f3(n,L), for a helium atom, which is 
based on another presentation of j3(n,L) in terms of the acceleration gauge dipole operators and a full diagonalization 
of the Hamiltonian. This method works well for two- and three-electron helium(lithium)-like atoms [4j but attempts 
O ■ to apply it to other systems like H^ molecular ion were not very successful. The other disadvantage of the method is 
a necessity to add some extra terms into a basis set, which looks like 1/r times regular solution. 

The major aim of a present work is to elaborate a universal method, which uses definition of the Bethe logarithm in 
terms of dipole operators in the velocity gauge and still as efficient as the Goldman-Drake method in case of hydrogen 
and helium like atoms. In fact, we have tried to carefully reanalyze the ideas of [2| and to present them in a more 
explicit and general form. The numerical results are confined to three-body systems while the theoretical expressions 
q ] are valid for a general few-body case. 

The paper is organized as follows. In Section II we consider a derivation of the leading order radiative corrections 
that allow us to define the Bethe logarithm for a general few-body Coulomb system as an integral over photon energy 

k. Next the asymptotic behaviour of the core integrand, J{k) — (j (Eg — H — h) J), and its first order perturbation 

wave function tpi(k), are derived, the leading order terms are obtained as expectation values of some operators. A 
variational property of the integral over k is discussed. It allows us to work out an efficient numerical scheme to 
calculate the Bethe logarithm using optimal parameters, which are variationally chosen. 

In Section III the numerical method, which is the main goal of our studies, is describe in details. And finally, the 
results of the Bethe logarithm calculation for the hydrogen, helium and Hj molecular ion ground states are presented 

^h ' and compared with other available results. 

t-H . We are using atomic units (h = e = m e = 1) throughout if something else is not explicitly stated. 

en : 

^ ! ii. theoretical overview 

^h ; 

^ . A. Radiative correction. Low energy contribution 

• i-H ' 

rS ■ Let us first consider a case of a hydrogen-like atom. 

The a(Za) 2 E nr order low-energy contribution, which results from the NRQED diagram (see Fig. [T](a)), may be 
written 
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where the last term is the "mass renormalization" contribution. Averaging the integrand over angular variables one 
gets 
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The integrand may be rewritten using the following operator identity 
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FIG. 1: NRQED diagrams for the self-energy and retarded transverse photon at ultra-low energies 



that results in 
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As was shown by Bcthc in 1947 [ff, the linearly divergent term should be associated with the "mass renormalization" 
of an electron and should be subtracted with the last term in expression ([3]). Thus, the remaining part may be splittcd 
onto a finite nonlogarithmic contribution 
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which results in appearance of the logarithmic term, the cut-of parameter is later canceled out by the logarithmic 
contribution from the high energy part. Here Eh is the Hartree energy. 

The high energy part is obtained from the one-loop scattering amplitude for an electron in an external field |6j 
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where E^" = (7^7^ — 7 l/ 7 M )/(2i), r = ln(A/m)+9/4— 21n(m/A mi „), and q 2 = 4m 2 sin 2 9. Here amplitude is expressed 
in the natural relativistic units (c = 1). 

At small g, scattering amplitude may be expanded: 
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The leading order contribution for a static scalar field with "renormalization" to a new infrared regularization 
parameter A, which is a cut-off of virtual quanta of momentum less than A (A = aA), is expressed 



r(0) 



a q 



Mi = -3^r2T + 6-8J ao+ 2^r- ■ 2 



a 1 



ta[q x pj 



CLQ 



(6) 



In order to get this expression, \nm/\ m i n should be replaced by [lnm/(2A) + 5/6] [6]. In the NRQED formalism this 
"renormalization" has been derived in [7|. 

In the coordinate space and atomic units the energy displacement due to the respective interaction is expressed in 
the NRQED by 
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Summing up the low energy contribution (|4|) and the high energy contribution ([7]), which comes from modification 
of the interaction of electron with a static field due to the one-loop self-energy diagram (Fig. [T](a)), one gets 
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where j3(n, I) is the so-called Bethe logarithm, it accumulates the contribution from (|4^) and will be formally defined 
later in Eq. (TTTT) . The cut-off parameter A cancels out. 

When both particles have finite masses one needs to include the retardation. The a(Za) 2 (m/M)E nr order low- 
energy contribution (Fig. [T](b)) may be written (here m a and mt are masses of the two particles) 
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The last term in the inner round brackets, — 1/fc, corresponds to the retardation interaction as it appears in the 
Breit-Pauli approximation (and is of order (Za) 2 (m/M)E nr ), and should be subtracted from the initial NRQED 
expression to avoid double counting. Consideration of the high-energy contribution (k > A), which comes from the 
same diagram, may be found in Ref. [8]. 

Summing up the contributions to the Bethe logarithm from Eqs. (]4ji) and ([9]) one finds that the dipole operator on 
the right and left-hand sides of the angle brackets may be recast as a nonrelativistic electric current density operator 
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as it may be expected, since in the quantum electrodynamics a virtual photon interacts with a current density. 

For a general many particle system the above speculations may be repeated directly in order to get the nonloga- 
rithmic part of the low-energy contribution. 

In summary, the Bethe logarithm may be defined as follows: 
Numerator: 
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And the Bethe logarithm itself is a ratio of these two quantities 



P(L,v) 



(11a) 



(lib) 



(lie) 



B. First order perturbation wave function, ipi(-), and asymptotic expansion of J(k) at k — > 00. 
The key quantity for our numerical studies is 

J{k)= (jiEo-H-ky 1 ^. (12) 
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Knowing this function one immediately gets a value for the nonrelativistic Bethe logarithm using Eq. (JTTJ). 
A general procedure to calculate J(k) is to solve an equation 

(E Q - H - k)tpi = iJipo, 



for different values of k. Since we are interested in asymptotic behaviour of J(k) for k — > 00, it is assumed that k is 
sufficiently large and as a first approximation one may take 
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Any approximate solution for ipi may be (formally) iteratively improved: 
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At small Tij, ipi should be smooth. In order to get a proper behaviour, one has to consider Eq. (TT3l) for r^ — > 
and keep only important terms 
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that gives homogeneous solutions of the type 



with fj>ij = y / 2rriijk. These solutions, taken for different pairs of particles, may be added to ij}\ to make the whole 
wave function smooth. So, we come to an approximation of ipi for k— >-oo, which is of required quality for our aims 
and has the following form, 



V4 



(1) 






i>3 



irii rru I r 



|. [1 _ e -mT« 



{1 + HijTij)] V'o(-)' 



(18) 



As is seen from this equation, there is no singular term in the wave function corresponding to a pair of identical 
particles. 

Integrand J(k) may be evaluated using the variational formalism as a stationary solution of a functional on t\j\ 

J(k) = -2 (V>o|a|V>i) - (MEo-H-k)^} . 

To get asymptotic expansion we substitute ip\ into this functional. A derivation of the asymptotic expansion for the 
hydrogen ground state and comparison with known analytical result may be found in Appendix A. 
At small r^ we get 
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where p_ means integration from to p. We assume that p satisfies l//ijj </)<!. 
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For regular r^ we use ip\ in a form: 



n 



(i) 



i 1 

T J ^0 + p [- H "^ J ]V'0 



Then 



and 



-2 (*|«j|tf>>- -f <j*>- | iMif^l 



i (HM) i ([if) . J]2> _ E 



fc 2 



fc 3 



p+ 



^ 2 m ij 






4n(8( rij )) 



That results in 
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Now we have to introduce a finite functional, which should replace a divergent expectation value of the 1/r 4 
operator: 
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where 
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Then summing up the Eqs. (fT§h ) and (119b ) one gets 
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For mixed terms: (ryrjy) / {rijTki) 2 in case of three-body calculation with the Hylleraas or exponential basis functions 
(see @, [l2[) a new type of singular integrals is required 
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A derivation of the explicit form for r_ 2i _ 2i o and stable recursions to get integrals for arbitrary n are presented in 
Appendix B. 



C. Variational property 



If we consider a quantity 
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we would find that for the ground state of a system this quantity possesses the variational property, since for the 
integrand for all k the following inequality is fulfilled 
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The same property remains satisfied for other states if integration is performed from some fco ~ 1: which lies above the 
poles related to the states E n < E$. It is known from the practical calculations that the low k contribution becomes 
numerically converged to a high accuracy at a moderate basis length of intermediate states, and thus with a good 
confidence the variational property, the higher the value of J7a the more accurate solution, is still remained in force. 
That allows us to perform optimization of the variational parameters of the basis set. 



III. NUMERICAL RESULTS 



A. Numerical scheme 



Here we consider the numerical scheme for the three-body Coulomb problem, which is then used in calculations of 
the Bethe logarithm for the helium and H 2 ground states. The wave functions both for the initial bound state and 



for the first order perturbation solution (or intermediate state), are taken in the form, 
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where 3^1(17, r 2 ) are the solid bipolar harmonics as defined in [llj . L is a total orbital angular momentum of a state. 
Complex parameters a, , /3j and 7$ are generated in a quasi-random manner jl2| | : 
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[a; J designates the fractional part of x, p a and q a are some prime numbers, [A\, A 2 ] and [A^A^] are real variational 
intervals which need to be optimized. Parameters fii and 7* are obtained in a similar way. 
Basis set for intermediate states is constructed as follows: 

1. First we use a regular basis set, which is taken similarly to the initial state with regular values of parameters 
(a, /?, 7) in exponentials. 

2. Then we build a special basis set with exponentially growing parameters for a particular r« 
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where r = A2/A1. 

Typically [A U A 2 ] = [2.5,4.5], and n max 



= 5 — 7, that corresponds to the photon energy interval k £ [0, 10 4 ]. 



3. For other pairs of (i, j) we take the similar basis sets as in 2. It is worthy to note that for identical particles 
this step should be omitted as is discussed in previous section. 

After the complete set of basis functions is constructed, we diagonalize matrix of the Hamiltonian Hi for intermediate 
states to get a set of (pseudo)state energies, E m , and then to calculate (0|iJ|to). These two sets of data are enough 
to restore J(fc): 
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From thus obtained J{k) we extrapolate coefficients of asymptotic expansion 
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which is taken in the same form as in analytic expression for the hydrogen atom (see Appendix A, Eq. (|A2jl ). The 
similar asymptotic expansion has been used in [1J]. The leading order terms of J(k) are obtained from Eq. (|21[) . 
That allows to get the high energy part of the numerator 
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N 



a 



p 



N 



40 2.98412855576549 

60 2.984128555765497 3 

80 2.98412855576549760 7 

100 2.9841285557654976107s 



2.9841284 9 45 

2.984128555i 66 

2.98412855575g 91 

2.9841285557654 4 120 



exact 



2.98412855576549761075977709002 



TABLE I: Convergence of the Bethe logarithm for the ground state of hydrogen and comparison with results of Drake and 
Goldman [jj]. 

Basis of the first type: 



[A_l, A_2] 
0.35000 0.80000 
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The total number of basis functions = 100 
Bethe logarithm ln(k_0) = 2.984128555765497610728 

FIG. 2: Parameters of the basis set and numerical value of the Bethe logarithm for the ground state of hydrogen for TV = 100. 

B. Results 



As a first example demonstrating capabilities of our method we take the ground state of a hydrogen atom. Results 
of numerical calculation are summarized in Table |TJ A basis set used for these calculations is a sum of exponentials with 
real parameters generated in a quasi-random way, the initial wave function is taken exact. On Fig. [2] an excerption 
of the listing of output with the parameters of the basis set for the case of N = 100 functions is shown. "Exact" in 
Table U means the value obtain by Huff's method jl4| . which is a series expansion. It is interesting to note that if the 
Huff series are taken with 100 terms only the result would be of the same accuracy as in our approach for N = 100. 
As is seen from this Table the new method demonstrates better accuracy than in Q, especially for small N, and good 
convergence rate. Worthy to note that variational parameters were roughly optimized for N = 80, and kept the same 
for all other cases. 

The main results of our studies are presented in Tables |H] and IIIII The first is the Bethe logarithm calculations for 
the ground state of helium and the next table is for the ground state of H^ molecular ion. Convergence is analyzed 
in two dimensions, with respect to the number of basis functions of the initial state and similarly for the basis size of 
the intermediate state. The major conclusion is that for ultimate precision it is important to check how the studied 
quantity depends on accuracy of the initial state. Indeed, it is reasonable to expect that the fractional error in the 
Bethe logarithm evaluation would be no better than the square root of the fractional uncertainty in the variational 
energy of the initial state, like in behaviour of expectation values of, say, 5-function or p 4 operators. 

Comparing our result for the ground state of helium with [HJ based on Goldman-Drake method (see Table IIV|) 
we see a discrepancy which requires to be somehow explained. This puzzle was a serious challenge for us, since the 
study of convergence (see Table II of [lfl]) unambiguously showed that the stability is achieved. So, we came back 
to our old calculations, which were based on the Goldman-Drake approach, and it was found that some essential 
part in the intermediate state wave function had been missed. If we add a new set of "regular" basis functions with 
exponentials having parameters a and f3 (in front of r\ and r<i, respectively) up to 40(!) then the numerical result 
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4.37016022306(2) 



TABLE II: Test of convergence of the Bethe logarithm value for the ground state of a helium atom. N a is a basis length for 
the initial state, Nb is a basis length for an intermediate state. 
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TABLE III: Test of convergence of the Bethe logarithm value for the ground state of a hydrogen molecular ion H J . N a is a 
basis length for the initial state, JV& is a basis length for an intermediate state. 

become j3 = 4.370 160 222 67 for the case of initial N — 1400 state used in [lj|. That showed a good coincidence with 
our "velocity gauge" result for the same initial state, as it might be expected! 

For convenience of comparison with other calculations we present here explicitly the asymptotic expansion of J{k) 
(k — > oo) for the helium ground state: 

1 /t2x _ 47r[(<S(ri)) + ((5(r 2 ))] fZ Z 2 ^2k - Z 3 In k + 20.00249948 

2 k (29) 

80.3063 Vk - 70.989 In k + 136.5 



J{k) i p2) = »Wl( 

k k z i 



Numerical coefficient in the first line has been calculated using Eg. (|21jl . Coefficients shown in a second line were 
obtained by a linear least squares fit using the SVD algorithm [l5|]. The fitting interval was k £ [20,6000], number of 
data points was 100, which were taken equidistant on the logarithmic scale of k, number of terms in the asymptotic 
expansion is 14-18. The results of fitting procedure are not sensitive to a number of data point. 
Similarly, the asymptotic expansion of J(k) for the ground state of a molecular ion H 2 is 
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where fa = m e Mi/ (m e +Mi) is the reduced mass of a respective proton of mass Mj. Since (<5(ri 2 )) < 10~ 10 , the term 
which is related to a nucleus-nucleus part (say, for HD + case) may be neglected. 

One additional remark on H 2 molecular ion is needed. Analyzing our new results for the Bethe logarithm we found 
that they have a systematic shift compared to the ones of [Ty] . That is because the formula for asymptotic expansion 
used in previous calculations (see Ref. [13], Eq. (12)) does not incorporate reduced masses into the leading terms of 
the expansion. This shift is almost state independent and (both for H^" [16| and HD + [13]) is of about 15 in the last 
two digits indicated in the Tables of [l6|, [I4] and should be subtracted. A new systematic calculation of the Bethe 
logarithm for the hydrogen molecular ions, which should provide better accuracy of about 8-9 significant digits, is in 
progress now. 
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Appendix A: Asymptotic expansion of J(k): Hydrogen ground state. 

Here the Hamiltonian and the charge current density operator are expressed 
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where ipo — 2Z 3 ' 2 e~ Zr is the ground state wave function and fi = v2fc, into the variational functional 
J(k) = -2 (V>o|V|Vi) ~ (MEo-H-k)^} 
• For r<p (p— ►(), and pp^>0) one gets: 
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Summing up, we obtain the leading terms of the asymptotic expansion: 
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k 



([H,V] 2 ) p+ -—4n(S(r)) 



Z 5 V2k Z 6 (lnfc-ln2 + l) 2Z 6 ( 7jB + lnp) 



k 3 k 3 k 3 

2Z 4 4Z 5 V2k 4Z 6 (lnfc-lnZ 2 ) 4Z 6 (31n2-l 

K* K**^ lC*^ K* 1 ^ 



Z- 3 4ir (6(r)) 



(Al) 
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which is exact up to free term in the 1/fc 3 order as it may be checked from comparison with analytical expression 
known for this case [19} (Z = l): 



J(k) = -384 



(1 + t)8(2-t) 



,fi(4 J 2-T,3-r;0 



1 2 4V2fc 41nfc-121n2+4 (21 + 2 7 r 2 )v / 2fc 81nfc-241n2+14+2r/>"(2) 



k k 2 k 3 



k 3 



3fc 4 



A: 4 



(A2) 



where r = Z/^/-2(E„ - k) and £ = [(1 - r)/(l + r)] 2 



Appendix B: r_ 2 ,-2,n(a,/3,7) 



In this appendix we show how to evaluate the integral 

r- 2 ,- 2 ,„(a,/3,7) = \ J J rr 2 r 2 - 2 ri l 2 e- OT1 -^-^ 2 d ri dr 2 dr 12 , 



(Bl) 



by means of 

/ d \ n 
r_ 2 ,_2,n(a, /9,7) = I ~ q- J r -2 ,-2,o (a, /3, 7)- 

To do this we have to obtain an analytical expression for r_2.-2,o(o^/3,7)- Let us consider 



r_ 2 _2,o(a,/3,7) 



1 



-ari 



drv 



dr 2 5" 

e2 ' 2 Jr!—r 2 



dr 12 e^ ri2 

e ~0 r 2 fr 1 +r 2 

dr 2 ~- 



rfr 



i2e 



I r°° e -(a+7)ri /•«> dr2 

= 7T" / « r l 9 / - S- I ' 



27J ei 



Ti '2 Jr 2 -ri 

-(/3-7)r 2 _ „-(/3+7)r 2 N 



/ dri 

27 4 



dri 



27J ei 



-(a+/9)n 



-(a+/3)n 



08 -7)- 



(0 + 7)- 



- (a+7)ri ffi((/3-7)n) 
rf 

"(°'-r)'-^i(C8+7)n) 



= 7r~ ^2(a+7, ei) L 2 (/3-7, e 2 ) - L 2 ((3+j, e 2 ) 
27 



+ — — I 2 [a + 7,/9-7, ei) — I 2 (a-7, P+7, d) 

27 27 



Here we use the notation from [2C 



-^-dt = a*- 1 r(l-p,a;c) ) 



I*(x, e) = V(l) - ln(xe), L 2 (x, e) = - - a; [^(2) - ln(ze)] - — , 



I P (x,y,e) 



t" 



■EMdt, 



I P (x,y,e) 



e xe Ei(ye) x 



-, ; ; In-i(x, v) L v (x + y, e), 

(p-l)eP" 1 p-l p u y> p-1 py y ' 



k 2 1 , 2 y 



1 2 

Ji(:c,2/,e) = -(lne+lny+7 B ) - — - -ln z - - dilog (1 



(B2) 



In expressions above i/j(ji) is a digamma function: ip(l) = — 7_e, ^>(n+l) = ip(n) + 1/n . 



11 



The integral is divergent and in order to get some meaningful finite expression one has to introduce some counterterm 
similar to what was done for the functional 1Z in Eq. (|20p. Still for our aims it is not needed, since we will be using 
T-2,-2,0 only for evaluation of finite integrals and, thus, any form of r_2.-2,o, which is self-consistent with other 
Fi mn (a, /3,7) via differentiation relations like (|B2I) . would be sufficient. So we may choose it in the following form: 



r_ 2> _ 2 ,o(a, P, 7) = ( 3 - 2 1E - — 1 7 - (a + p) 

a + l / n \\ 2 P + l , 



2 :(2- 7B -ln(a+7)) - ^(2 - 1e ~ HP+l)) 
+(a + p) ln(a + 0) - (a + 7) ln(a + 7) - (0 + 7) ln(/3 + 7) 

a/37r 2 



(B3) 



(a + 7)(/3+7) ln 2 (<* + \ 



4", 

(a- 7 )(/3 + 7 ) 
27 



dilog 



P + l 
a + [3 

J+l 



12 7 

Q9-7)(tt + 7) 

2 7 



dilog 



/? 



a + 7 



Evaluation of r_ 2 ,-2.n is straightforward (see [2l| for details) except for two terms, which require additional remarks. 
We introduce two recursions: 



An = -■=- 



j9 

c*- 



and 



£n = -7T- 



^7 



— dilog 

7 



7 



•0 



a + 7 



a + 7 



- [nAn-x 

7 



7 



Bn] 



Bn = — =- 



j9 



dilog 



/3 



a + 7 



F„ = - — 



d-, 



hr 



a + 7 

J+l 



(B4) 



(B5) 



At first glance these recursions are not stable, when 7 is small. However, more careful analysis shows that, say, A n —\ 
and B n (as well as E n -i and F n ) are of the same sign, and no subtraction, which leads to loss of numerical accuracy, 
occurs. 
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